/*

PhyML:  a program that  computes maximum likelihood phyLOGenies from
DNA or AA homoLOGous sequences.

Copyright (C) Stephane Guindon. Oct 2003 onward.

All parts of the source except where indicated are distributed under
the GNU public licence. See http://www.opensource.org for details.

*/

#include "cl.h"


/*********************************************************/
/**
* Fill the Option fields, with the argc array
*/
void Read_Command_Line(option *io, int argc, char **argv)
{
  int c;
  int open_ps_file;
  int use_gamma;
  int writemode;

  if(argc == 1) Exit("\n. No argument was passed to the program. Please check the documentation. \n");

  struct option longopts[] =
    {
      {"n_rgrft",           required_argument,NULL,0},
      {"n_globl",           required_argument,NULL,1},
      {"max_dist",          required_argument,NULL,2},
      {"n_optim",           required_argument,NULL,3},
      {"n_best",            required_argument,NULL,4},
      {"model",             required_argument,NULL,5},
      {"search",            required_argument,NULL,6},
      {"datatype",          required_argument,NULL,7},
      {"multiple",          required_argument,NULL,8},
      {"input",             required_argument,NULL,9},
      {"bootstrap",         required_argument,NULL,10},
#ifdef ENABLE_EMPIRICAL_BAYES
      {"support",           required_argument,NULL,10}, // same as bootstrap
#endif
      {"ts/tv",             required_argument,NULL,11},
      {"nclasses",          required_argument,NULL,12},
      {"pinv",              required_argument,NULL,13},
      {"alpha",             required_argument,NULL,14},
      {"inputtree",         required_argument,NULL,15},
      {"min_diff_lk_local", required_argument,NULL,16},
      {"min_diff_lk_global",required_argument,NULL,17},
      {"steph_spr",         no_argument,NULL,18},
      {"brent_it_max",      required_argument,NULL,19},
      {"rand_start",        no_argument,NULL,20},
      {"n_rand_starts",     required_argument,NULL,21},
      {"sequential",        no_argument,NULL,22},
      {"inside_opt",        no_argument,NULL,23},
      {"p_moves",           required_argument,NULL,24},
      {"fast_nni",          no_argument,NULL,25},
      {"g_pars",            no_argument,NULL,26},
      {"r_seed",            required_argument,NULL,27},
      {"collapse_boot",     required_argument,NULL,28},
      {"random_boot",       required_argument,NULL,29},
      {"print_trace",       no_argument,NULL,30},
      {"print_site_lnl",    no_argument,NULL,31},
      {"cov",               no_argument,NULL,32},
      {"cov_delta",         required_argument,NULL,33},
      {"cov_alpha",         required_argument,NULL,34},
      {"cov_ncats",         required_argument,NULL,35},
      {"ps",                no_argument,NULL,36},
      {"cov_free",          no_argument,NULL,37},
      {"no_gap",            no_argument,NULL,38},
      {"n_rr_branch",       required_argument,NULL,39},
      {"append",            no_argument,NULL,40},
      {"no_five_branch",    no_argument,NULL,41},
      {"pars_thresh",       required_argument,NULL,42},
      {"min_diff_lk_move",  required_argument,NULL,43},
      {"hybrid",            no_argument,NULL,44},
      {"use_median",        no_argument,NULL,45},
      {"run_id",            required_argument,NULL,46},
      {"pars",              no_argument,NULL,47},
      {"quiet",             no_argument,NULL,48},
      {"version",           no_argument,NULL,49},
      {"calibration",       required_argument,NULL,50},
      {"boot_progress_every", required_argument,NULL,51},
      {"aa_rate_file",        required_argument,NULL,52},
      {"chain_len",           required_argument,NULL,53},
      {"sample_freq",         required_argument,NULL,54},
      {"burnin",              required_argument,NULL,55},
      {"no_memory_check",     no_argument,NULL,56},
      {"no_colalias",         no_argument,NULL,57},
      {"alias_subpatt",       no_argument,NULL,58},
      {"no_data",             no_argument,NULL,59},
      {"fastlk",              required_argument,NULL,60},
#ifdef ENABLE_MULTIMAX
      {"opt_algorithm",       required_argument,NULL,1061},
      {"bfgs_stepsize",       required_argument,NULL,1062},
      {"bfgs_itmax",          required_argument,NULL,1063},
      {"lnsearch_stepmax",    required_argument,NULL,1064},
      {"swap_before_opt",	  required_argument,NULL,1065},
      {"nt_rates",	          required_argument,NULL,1066},
      {"print_ml_delta", 	  required_argument,NULL,1067},
      {"print_nni_scores", 	  required_argument,NULL,1068},
      {"opt_subst_param",     required_argument,NULL,1069},
      {"opt_topo",     		  required_argument,NULL,1071},
      {"first_branch",     		required_argument,NULL,1072},
      {"l_max",     		required_argument,NULL,1073},
#endif
#ifdef ENABLE_EMPIRICAL_BAYES
      {"eb_ngens",            required_argument,NULL,1090},
#endif
#ifdef ENABLE_FULL_SWAP_OPT
      {"full_swap_opt",     required_argument,NULL,1070},
#endif
#ifdef ENABLE_MPI_TREES
      {"parallelism_type",  required_argument,NULL,1080},
#endif
#ifdef ENABLE_MULTIMAX
		{"sa_iters_per_stage",   required_argument,NULL,2000},
		{"sa_num_anneal_stages", required_argument,NULL,2001},
		{"sa_scale_temp",         required_argument,NULL,2002},
		{"sa_temp_end",          required_argument,NULL,2003},
		{"sa_set_back",          required_argument,NULL,2004},
		{"sa_temp_start",        required_argument,NULL,2005},
		{"sa_max_brlen",       required_argument,NULL,2006},
		{"sa_max_pinvar",      required_argument,NULL,2007},
		{"sa_max_gamma",       required_argument,NULL,2008},
		{"sa_max_kappa",       required_argument,NULL,2009},
		{"sa_max_lambda",       required_argument,NULL,2010},
		{"sa_max_blprop",			required_argument,NULL,2011},
		{"sa_max_pi",			required_argument,NULL,2012},
		{"sa_max_rr",			required_argument,NULL,2013},
		{"sa_max_emig",        required_argument,NULL,2014},
		{"sa_prob_nni",          required_argument,NULL,2015},
		{"sa_prob_spr",          required_argument,NULL,2016},
		{"sa_prob_brlen",        required_argument,NULL,2017},
		{"sa_prob_gamma",        required_argument,NULL,2018},
		{"sa_prob_kappa",        required_argument,NULL,2019},
		{"sa_prob_lambda",       required_argument,NULL,2020},
		{"sa_prob_rr",           required_argument,NULL,2021},
		{"sa_prob_blprop",		required_argument,NULL,2022},
		{"sa_prob_topology",     required_argument,NULL,2023},
		{"sa_prob_pinvar",       required_argument,NULL,2024},
		{"sa_prob_pi",           required_argument,NULL,2025},
		{"sa_prob_emig",         required_argument,NULL,2026},
		{"sa_tau_start",		  required_argument,NULL,2027},
		{"sa_tau_end",			  required_argument,NULL,2028},
		{"sa_sele_bl",		  	  required_argument,NULL,2029},
		{"sa_sele_blprop",		  required_argument,NULL,2030},
		{"sa_sele_gamma",		  required_argument,NULL,2031},
		{"sa_sele_pinvar",		  required_argument,NULL,2032},
		{"sa_sele_kappa",		  required_argument,NULL,2033},
		{"sa_sele_lambda",		  required_argument,NULL,2034},
		{"sa_sele_rr",		  	  required_argument,NULL,2035},
		{"sa_sele_pi",		  	  required_argument,NULL,2036},
		{"sa_preopt",			  required_argument,NULL,2037},
		{"sa_sele_max_params",	  required_argument,NULL,2038},
		{"sa_stop_early",         required_argument,NULL,2039},

		// These are for drawing the local likelihod landscape of branch lengths
		{"debug_show_walk",			required_argument,NULL,2050},
		{"debug_n_steps_walk",		required_argument,NULL,2051},
		{"stepsize_walk",           required_argument,NULL,2052},

		{"print_bl_trace",           required_argument,NULL,2053},

		// For mapping likelihood landscapes in my publication, August 2011.
		{"fix_param",				required_argument,NULL,2060},
		{"fix_param_min_val",		required_argument,NULL,2061},
		{"fix_param_max_val",		required_argument,NULL,2062},
		{"fix_param_step",			required_argument,NULL,2063},
#endif
#ifdef ENABLE_MBL
		{"bl_mixtures",         	required_argument,NULL,3000}, // the number of BL mixtures
		{"bl_mixture_wts",         required_argument,NULL,3001}, // the weights on those mixtures
		{"opt_bl_mixtures",        required_argument,NULL,3002}, // should we optimize those values?
#endif
#ifdef ENABLE_ASR
		{"asr",        required_argument,NULL,3030},
#endif
      {0,0,0,0}
    };


  writemode = 1;
  open_ps_file = 0;
  use_gamma = 0;
  while((c = getopt_long(argc,argv,"qi:d:m:b:n:t:f:zk:v:c:a:u:ho:s:x:g:l:ep",longopts,NULL)) != -1)
    {
      switch(c)
	{

#ifdef ENABLE_EMPIRICAL_BAYES
      case 1090:
      {
		int tmp = io->eb_ngens;
		io->eb_ngens = atoi(optarg);
		if(io->eb_ngens < 1) io->eb_ngens = tmp;
		break;
      }
#endif

#ifdef ENABLE_MULTIMAX
		case 2000:
		{
			int tmp = io->sa_iters_per_stage;
			io->sa_iters_per_stage = atoi(optarg);
			if(io->sa_iters_per_stage < 1) io->sa_iters_per_stage = tmp;
			break;
		}
		case 2001:
		{
			int tmp = io->sa_num_anneal_stages;
			io->sa_num_anneal_stages = atoi(optarg);
			if(io->sa_num_anneal_stages < 1) io->sa_num_anneal_stages = tmp;
			break;
		}
		case 2002:
		{
			int tmp = io->sa_scale_temp;
			io->sa_scale_temp = atoi(optarg);
			if(io->sa_scale_temp != 0 && io->sa_scale_temp != 1) io->sa_scale_temp = tmp;
			break;
		}
		case 2003:
		{
			double tmp = io->sa_temp_end;
			io->sa_temp_end = atof(optarg);
			if(io->sa_temp_end < 0.0) io->sa_temp_end = tmp; // use the default T, if the user's argument was < 0.0
			if(io->sa_temp_end == 0.0) io->sa_temp_end = 0.00000001; // the final T cannot actually be 0.0, otherwise we get error when calcuting the Boltzmann's P.
			break;
		}
		//		{"sa_temp_end",          required_argument,NULL,55},
		case 2004:
		{
			int tmp = io->sa_set_back;
			io->sa_set_back = atoi(optarg);
			if(io->sa_set_back < 0) io->sa_set_back = tmp;
			break;
		}
		//		{"sa_set_back",          required_argument,NULL,56},
		case 2005:
		{
			double tmp = io->sa_temp_start;
			io->sa_temp_start = atof(optarg);
			if(io->sa_temp_start < 0.0) io->sa_temp_start = tmp;
			break;
		}
		case 2006:
		{
			double tmp = io->sa_sigma_bl;
			io->sa_sigma_bl = atof(optarg);
			if(io->sa_sigma_bl < 0.0) io->sa_sigma_bl = tmp;
			break;
		}
		//		{"sa_sigma_bl",       required_argument,NULL,59},
		case 2007:
		{
			double tmp = io->sa_sigma_pinvar;
			io->sa_sigma_pinvar = atof(optarg);
			if(io->sa_sigma_pinvar < 0.0) io->sa_sigma_pinvar = tmp;
			break;
		}
		//		{"sa_sigma_pinvar",      required_argument,NULL,60},
		case 2008:
		{
			double tmp = io->sa_sigma_gamma;
			io->sa_sigma_gamma = atof(optarg);
			if(io->sa_sigma_gamma < 0.0) io->sa_sigma_gamma = tmp;
			break;

		}
		//		{"sa_sigma_gamma",       required_argument,NULL,61},
		case 2014:
		{
			double tmp = io->sa_sigma_emig;
			io->sa_sigma_emig = atof(optarg);
			if(io->sa_sigma_emig < 0.0) io->sa_sigma_emig = tmp;
			break;
		}
		case 2009:
		{
			double tmp = io->sa_sigma_kappa;
			io->sa_sigma_kappa = atof(optarg);
			if(io->sa_sigma_kappa < 0.0) io->sa_sigma_kappa = tmp;
			break;
		}
		case 2010:
		{
			double tmp = io->sa_sigma_lambda;
			io->sa_sigma_lambda = atof(optarg);
			if(io->sa_sigma_lambda < 0.0) io->sa_sigma_lambda = tmp;
			break;
		}
		case 2011:
		{
			double tmp = io->sa_sigma_blprop;
			io->sa_sigma_blprop = atof(optarg);
			if(io->sa_sigma_blprop < 0.0) io->sa_sigma_blprop = tmp;
			break;
		}
		case 2012:
		{
			double tmp = io->sa_sigma_pi;
			io->sa_sigma_pi = atof(optarg);
			if(io->sa_sigma_pi < 0.0) io->sa_sigma_pi = tmp;
			break;
		}
		case 2013:
		{
			double tmp = io->sa_sigma_rr;
			io->sa_sigma_rr = atof(optarg);
			if(io->sa_sigma_rr < 0.0) io->sa_sigma_rr = tmp;
			break;
		}
		case 2037:
		{
			// -1 = no pre-optimization
			// 0 = pre-opt using SLM
			// 1 = pre-opt using BFGS
			int tmp = io->sa_preopt;
			io->sa_preopt = atoi(optarg);
			if(io->sa_preopt < 0 && io->sa_preopt > 2) io->sa_preopt = tmp;
			break;
		}
		case 2039:
		{
		  int tmp = io->sa_stop_early;
		  io->sa_stop_early = atoi(optarg);
		  if(io->sa_stop_early < 0 && io->sa_stop_early > 2) io->sa_stop_early = tmp;
		  break;
		}

		case 2050:
		{
			  int tmp = io->show_walk;
			  io->show_walk = atoi(optarg);
			  if(io->show_walk < 0 && io->show_walk > 3) io->show_walk = tmp;
			  break;
		}
		case 2051:
		{
			  int tmp = io->n_steps_walk;
			  io->n_steps_walk = atoi(optarg);
			  if(io->n_steps_walk < 0) io->n_steps_walk = tmp;
			  break;
		}
		case 2052:
		{
			  phydbl tmp = io->stepsize_walk;
			  io->stepsize_walk = atof(optarg);
			  if(io->stepsize_walk < 0.0) io->stepsize_walk = tmp;
			  break;
		}
		case 2053:
		{
			  phydbl tmp = io->print_bl_trace;
			  io->print_bl_trace = atof(optarg);
			  if(io->print_bl_trace < 0.0) io->print_bl_trace = tmp;
			  break;
		}
		case 2060:
		{
			  int tmp = io->fix_param;
			  io->fix_param = atoi(optarg);
			  if(io->fix_param < 0) io->fix_param = tmp;
			  break;
		}
		case 2061:
		{
			  phydbl tmp = io->fix_param_min_val;
			  io->fix_param_min_val = atof(optarg);
			  if(io->fix_param_min_val < 0.0) io->fix_param_min_val = tmp;
			  break;
		}
		case 2062:
		{
			  phydbl tmp = io->fix_param_max_val;
			  io->fix_param_max_val = atof(optarg);
			  if(io->fix_param_max_val < 0.0) io->fix_param_max_val = tmp;
			  break;
		}
		case 2063:
		{
			  phydbl tmp = io->fix_param_step;
			  io->fix_param_step = atof(optarg);
			  if(io->fix_param_step < 0.0) io->fix_param_step = tmp;
			  break;
		}

		case 2015:
		{
			double tmp = io->sa_prob_NNI;
			io->sa_prob_NNI = atof(optarg);
			if(io->sa_prob_NNI < 0.0) io->sa_prob_NNI = tmp;
			break;
		}
		//		{"prob_nni",          required_argument,NULL,63},
		case 2016:
		{
			double tmp = io->sa_prob_SPR;
			io->sa_prob_SPR = atof(optarg);
			if(io->sa_prob_SPR < 0.0) io->sa_prob_SPR = tmp;
			break;
		}
		//		{"prob_spr",          required_argument,NULL,64},
		case 2017:
		{
			double tmp = io->sa_prob_brlen;
			io->sa_prob_brlen = atof(optarg);
			if(io->sa_prob_brlen < 0.0) io->sa_prob_brlen = tmp;
			break;
		}
		//		{"sa_prob_brlen",        required_argument,NULL,65},
		case 2018:
		{
			double tmp = io->sa_prob_gamma;
			io->sa_prob_gamma = atof(optarg);
			if(io->sa_prob_gamma < 0.0) io->sa_prob_gamma = tmp;
			break;
		}
		//		{"sa_prob_gamma",        required_argument,NULL,66},
		case 2019:
		{
			double tmp = io->sa_prob_kappa;
			io->sa_prob_kappa = atof(optarg);
			if(io->sa_prob_kappa < 0.0) io->sa_prob_kappa = tmp;
			break;

		}
		//		{"sa_prob_kappa",        required_argument,NULL,67},
		case 2020:
		{
			double tmp = io->sa_prob_lambda;
			io->sa_prob_lambda = atof(optarg);
			if(io->sa_prob_lambda < 0.0) io->sa_prob_lambda = tmp;
			break;
		}
		//		{"sa_prob_lambda",       required_argument,NULL,68},
		case 2021:
		{
			double tmp = io->sa_prob_rr;
			io->sa_prob_rr = atof(optarg);
			if(io->sa_prob_rr < 0.0) io->sa_prob_rr = tmp;
			break;
		}
		//		{"sa_prob_rr",           required_argument,NULL,69},
		case 2022:
		{
			double tmp = io->sa_prob_blprop;
			io->sa_prob_blprop = atof(optarg);
			if(io->sa_prob_blprop < 0.0) io->sa_prob_blprop = tmp;
			break;
		}

		case 2029:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_bl = x;
			}
			break;
		}
		case 2030:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1 || x == 2){
				io->sa_sele_blprop = x;
			}
			break;
		}
		case 2031:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_gamma = x;
			}
			break;
		}
		case 2032:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_pinvar = x;
			}
			break;
		}
		case 2033:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_kappa = x;
			}
			break;
		}
		case 2034:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_lambda = x;
			}
			break;
		}
		case 2035:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1 || x == 2){
				io->sa_sele_rr = x;
			}
			break;
		}
		case 2036:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1 || x == 2){
				io->sa_sele_pi = x;
			}
			break;
		}
		case 2023:
		{
			double tmp = io->sa_prob_topo;
			io->sa_prob_topo = atof(optarg);
			if(io->sa_prob_topo < 0.0) io->sa_prob_topo = tmp;
			break;
		}
		//		{"sa_prob_topo",     required_argument,NULL,71},
		case 2024:
		{
			double tmp = io->sa_prob_pinvar;
			io->sa_prob_pinvar = atof(optarg);
			if(io->sa_prob_pinvar < 0.0) io->sa_prob_pinvar = tmp;
			break;
		}
		//		{"sa_prob_pinvar",       required_argument,NULL,72},
		case 2025:
		{
			double tmp = io->sa_prob_pi;
			io->sa_prob_pi = atof(optarg);
			if(io->sa_prob_pi < 0.0) io->sa_prob_pi = tmp;
			break;
		}
		case 2038:
		{
			int tmp = io->sa_sele_max_params;
			int x = atoi(optarg);
			if (x > 0){
				io->sa_sele_max_params = x;
			}
			else
				io->sa_sele_max_params = tmp;
			break;
		}
#endif


#ifdef ENABLE_MPI_TREES
      case 1080:
      {
    	  int tmp = io->parallelism_type;
    	  io->parallelism_type = atoi(optarg);
    	  if (io->parallelism_type < 0 || io->parallelism_type > 2) io->parallelism_type = tmp;
    	  break;
      }
#endif

#ifdef ENABLE_FULL_SWAP_OPT
	  case 1070:
	  {
			int tmp = io->full_swap_opt;
			io->full_swap_opt = atoi(optarg);
			if (io->full_swap_opt < 0 || io->full_swap_opt > 1) io->full_swap_opt = tmp;
			break;
	  }
#endif
#ifdef ENABLE_MULTIMAX
	  case 1065:
	  {
			int tmp = io->swap_before_opt;
			io->swap_before_opt = atoi(optarg);
			if(io->swap_before_opt < 0) io->swap_before_opt = tmp;
			break;
	  }
	  case 1064:
	  {
			double tmp = io->lnsearch_stepmax;
			io->lnsearch_stepmax = atof(optarg);
			if (io->lnsearch_stepmax < 0.0) io->lnsearch_stepmax = tmp;
			break;
	  }
      case 1063:
      {
			int tmp = io->bfgs_itmax;
			io->bfgs_itmax = atoi(optarg);
			if(io->bfgs_itmax < 0) io->bfgs_itmax = tmp;
			break;
      }
	  case 1062:
	  {
			double tmp = io->bfgs_stepsize;
			io->bfgs_stepsize = atof(optarg);
			if (io->bfgs_stepsize < 0.0) io->bfgs_stepsize = tmp;
			break;
	  }
      case 1061:
      {
			int tmp = io->opt_algorithm;
			io->opt_algorithm = atoi(optarg);
			if(io->opt_algorithm < 0) io->opt_algorithm = tmp;
			break;
      }
#endif
#ifdef ENABLE_MBL
      case 3000: // specify the number of BL mixture classes
      {
              int user_n_l = atoi(optarg);
              if ( (n_l > 1) && (user_n_l != n_l) ) // if the user previously specified the size of the mixture model (probably using --bl_mixture_wts)...
              {
				  PhyML_Printf ("\n. Yikes! I found %d mixture weights after the option --bl_mixture_wts, but you specified %d mixtures after the option --bl_mixtures.  I am confused!", n_l, user_n_l);
				  PhyML_Printf("\n. Type any key to exit.\n");
				  Exit("\n");
              }
              n_l = user_n_l;

              if(n_l < 1){
                      PhyML_Printf("\n. The number of branch length mixture classes (--bl_mixtures) must be at least 1.");
                      PhyML_Printf("\n. Type any key to exit.\n");
                      Exit("\n");
              }
              break;
      }
      case 3001: // BL mixture weights
	  {
			 // io->mod->s_opt->user_bl_wts = 1; // yes, we're using user-specified BL mixture weights

			  // 1. Sanity check:
			  if(optarg[0] != '['){
					  char choix;
					  PhyML_Printf ("\n. The branch length mixture weights must be in the form [NUM1,NUM2,...,NUMN] surrounded by brackets '[' and ']', separated by commas, and with not spaces");
					  PhyML_Printf("\n. Type any key to exit.\n");
					  if(!scanf("%c",&choix)) Exit("\n");
					  Exit("\n");
			  }

			  // 2. determine how many BL categories are implied
			  // Although its tempting to just query the global n_l, that parameter
			  // is not guaranteed to have been parsed from the CLI yet.
			  int count_mixtures = 0; // use this as the counter for the number of props read in.
			  int i = 0;
			  for(i=1; i < strlen(optarg); i++)
			  {
					  if (optarg[i] == ',')
					  {
							  count_mixtures++;
					  }
					  if (optarg[i] == ']')
					  {
							  count_mixtures++;
							  break;
					  }
			  }

			  // 2b. Sanity check:
			  if (count_mixtures != n_l)
			  {
				  if (n_l > 1) // if the user specified a discordant number of mixtures...
				  {
					  PhyML_Printf ("\n. Yikes! I found %d mixture weights after the option --bl_mixture_wts, but you specified %d mixtures after the option --bl_mixtures.  I am confused!", count_mixtures, n_l);
					  PhyML_Printf("\n. Type any key to exit.\n");
					  Exit("\n");
				  }
				  else // else the user did not specify any number of mixtures...
				  {
					  PhyML_Printf("\n. --> I detected %d mixture weights after the option --bl_mixture_wts.  I will now enable the mixture model to have %d mixtures.\n", count_mixtures, count_mixtures);
					  n_l = count_mixtures;
				  }
			  }

			  // 3. Allocate memory for 'count_mixtures' number of BL sets
			  io->mod->bl_mixture_wts = (phydbl *)mCalloc(count_mixtures, sizeof(phydbl));

              // 4. finally, parse the BL proportions as floating-point values...
			  count_mixtures = 0; // reset counter
              i = 0;           // reset counter
              int j = 0;
              char aprop[100];
              for(i=1; i < strlen(optarg); i++)
              {
                      if(optarg[i] != ',' && optarg[i] != ']')
                      {
                              aprop[j] = optarg[i];
                              j++;
                      }
                      else
                      {
                              aprop[j+1] = '\0'; // end the string
                              j = 0; //reset j
                              io->mod->bl_mixture_wts[count_mixtures] = atof(aprop);
                              count_mixtures++;
                              //PhyML_Printf("\n. cl.c 621 mixture %d wt = %f", count_mixtures, atof(aprop) );
                      }
              }
              Set_nth_mixture_weight(io->mod->bl_mixture_wts, count_mixtures);
              break;
	  }
      case 3002: // should we optimize the BL mixture weights?
      {
			int tmp = io->mod->s_opt->opt_bl_mixtures;
			io->mod->s_opt->opt_bl_mixtures = atoi(optarg);
			if((io->mod->s_opt->opt_bl_mixtures < 0) || (io->mod->s_opt->opt_bl_mixtures > 1)) io->mod->s_opt->opt_bl_mixtures = tmp;
			break;
      }
#endif

#ifdef ENABLE_ASR
      case 3030:
	  {
		  int tmp = io->asr_algo;
		  io->asr_algo = atoi(optarg);
		  if(io->asr_algo < 0 || io->asr_algo > 2) io->asr_algo = tmp;
		  break;
	  }
#endif
	case 60:
	  {
	    if(!strcmp(optarg,"yes") || !strcmp(optarg,"YES"))
	      {
		io->lk_approx = NORMAL;
	      }
	    else if(!strcmp(optarg,"no") || !strcmp(optarg,"NO"))
	      {
		io->lk_approx = EXACT;
	      }
	    else
	      {
		PhyML_Printf("\n. 'fastlk' option: please use 'yes' or 'no'");
		Exit("\n");
	      }

	    break;
	  }
	case 59:
	  {
	    io->mcmc->use_data = NO;
	    break;
	  }
	case 58:
	  {
	    io->do_alias_subpatt = YES;
	    break;
	  }
	case 57:
	  {
	    io->colalias = NO;
	    break;
	  }
	case 56:
	  {
	    io->mem_question = NO;
	    break;
	  }
	case 55:
	  {
	    phydbl len;
	    len = atof(optarg);
	    io->mcmc->chain_len_burnin = (int)len;
	    if(io->mcmc->chain_len_burnin < 1)
	      {
		char choix;
		PhyML_Printf("\n. chain_len_burnin must be an integer greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case 54:
	  {
	    phydbl len;
	    len = atof(optarg);
	    io->mcmc->sample_interval = (int)len;
	    if(io->mcmc->sample_interval < 1)
	      {
		char choix;
		PhyML_Printf("\n. sample_interval must be an integer greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case 53:
	  {
	    phydbl len;
	    len = atof(optarg);
	    io->mcmc->chain_len = (int)len;
	    if(io->mcmc->chain_len < 1)
	      {
		char choix;
		PhyML_Printf("\n. chain_len must be an integer greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case 52:
	  {
	    char *s;
	    s = (char *)mCalloc(T_MAX_FILE, sizeof(char));
	    strcpy(s,optarg);
	    io->fp_aa_rate_mat = Openfile(s,0);
	    Free(s);
	    break;
	  }
	case 51:
	  {
	    io->boot_prog_every = atoi(optarg);
	    if(io->boot_prog_every < 1)
	      {
		char choix;
		PhyML_Printf("\n. boot_progress_every must be an integer greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case 50:
	  {
	    strcpy(io->clade_list_file,optarg);
	    break;
	  }
	case 49:
	  {
	    PhyML_Printf("%s\n",VERSION);
	    Exit("");
	    break;
	  }
	case 48 :
	  {
	    io->quiet = 1;
	    break;
	  }
	case 'p' : case 47 :
	  {
	    io->in_tree = 1;
	    break;
	  }
	case 46 :
	  {
	    io->appebr_run_ID = 1;
	    strcpy(io->run_id_string,optarg);
	    break;
	  }
	case 45 :
	  {
	    io->mod->gamma_median = 1;
	    break;
	  }
	case 44 :
	  {
	    io->mod->s_opt->hybrid_thresh = 0;
	    break;
	  }
	case 43 :
	  {
	    io->mod->s_opt->min_diff_lk_move = atof(optarg);
	    if(io->mod->s_opt->min_diff_lk_move < 0)
	      {
		char choix;
		PhyML_Printf("\n. Min_diff_lk_move must be a double greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case 42 :
	  {
	    io->mod->s_opt->pars_thresh = (int)atoi(optarg);
	    if(io->mod->s_opt->pars_thresh < 0)
	      {
		PhyML_Printf("\n. The parsimony threshold must be an integer greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		Exit("\n");
	      }
	    break;
	  }
	case 41 :
	  {
	    io->mod->s_opt->opt_five_branch = 0;
	    break;
	  }
	case 40 :
	  {
	    writemode = 2;
	    break;
	  }
	case 39 :
	  {
	    char choix;
	    io->mod->n_rr_branch = (int)atoi(optarg);
	    if(io->mod->n_rr_branch < 1)
	      {
		PhyML_Printf("\n. The number of classes must be an integer greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case 38 :
	  {
	    io->rm_ambigu = 1;
	    break;
	  }
	case 37 :
	  {
	    io->mod->s_opt->opt_cov_free_rates = 1;
	    io->mod->m4mod->use_cov_alpha      = 0;
	    io->mod->m4mod->use_cov_free       = 1;
	    break;
	  }
	case 36 :
	  {
	    open_ps_file = 1;
	    break;
	  }
	case 35 :
	  {
	    io->mod->m4mod->n_h = (int)atoi(optarg);

	    if(io->mod->m4mod->n_h < 1)
	      {
		char choix;
		PhyML_Printf("\n. The number of classes must be greater than 0.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case 34 :
	  {
	    io->mod->m4mod->use_cov_alpha = 1;
	    io->mod->m4mod->use_cov_free  = 0;

	    if(!strcmp(optarg,"e") || !strcmp(optarg,"E") ||
	       !strcmp(optarg,"estimated") || !strcmp(optarg,"ESTIMATED"))
	      {
		io->mod->s_opt->opt_cov_alpha = 1;
		io->mod->m4mod->alpha         = 1.0;
	      }
	    else
	      {
		io->mod->m4mod->alpha = (phydbl)atof(optarg);

		if(io->mod->m4mod->alpha < 1.E-5)
		  {
		    char choix;
		    PhyML_Printf("\n. The value of alpha must be greater than 1.E-5.\n");
		    PhyML_Printf("\n. Type any key to exit.\n");
		    if(!scanf("%c",&choix)) Exit("\n");
		    Exit("\n");
		  }
	      }
	    break;
	  }
	case 33 :
	  {
	    if(!strcmp(optarg,"e") || !strcmp(optarg,"E") ||
	       !strcmp(optarg,"estimated") || !strcmp(optarg,"ESTIMATED"))
	      {
		io->mod->s_opt->opt_cov_delta = 1;
		io->mod->m4mod->delta         = 1.0;
	      }
	    else
	      {
		io->mod->m4mod->delta = (phydbl)atof(optarg);

		if(atof(optarg) < 1.E-10)
		  {
		    char choix;
		    PhyML_Printf("\n. The value of delta must be larger than 1.E-10.\n");
		    PhyML_Printf("\n. Type any key to exit.\n");
		    if(!scanf("%c",&choix)) Exit("\n");
		    Exit("\n");
		  }
	      }
	    break;
	  }
	case 32 :
	  {
	    io->mod->use_m4mod = YES;
	    break;
	  }
	case 31 :
	  {
	    io->print_site_lnl = 1;
	    break;
	  }
	case 30 :
	  {
	    io->print_trace = 1;
	    break;
	  }
	case 29 :
	  {
	    io->random_boot_seq_order = (int)atoi(optarg);
	    break;
	  }
	case 28 :
	  {
	    io->collapse_boot = (int)atoi(optarg);
	    break;
	  }
	case 27 :
	  {
	    io->r_seed = (int)atoi(optarg);
	    break;
	  }
	case 26 :
	  {
	    io->mod->s_opt->general_pars = 1;
	    break;
	  }
	case 25 :
	  {
	    io->mod->s_opt->fast_nni = 1;
	    break;
	  }
	case 24 :
	  {
	    io->mod->s_opt->p_moves_to_examine = (phydbl)atof(optarg);
	    break;
	  }
	case 23 :
	  {
	    io->mod->s_opt->wim_inside_opt = 1;
	    break;
	  }
	case 0 :
	  {
	    io->mod->s_opt->wim_n_rgrft = atoi(optarg);
	    break;
	  }
	case 1 :
	  {
	    io->mod->s_opt->wim_n_globl = atoi(optarg);
	    break;
	  }
	case 2 :
	  {
	    io->mod->s_opt->wim_max_dist = atoi(optarg);
	    break;
	  }
	case 3 :
	  {
	    io->mod->s_opt->wim_n_optim = atoi(optarg);
	    break;
	  }
	case 4 :
	  {
	    io->mod->s_opt->wim_n_best = atoi(optarg);
	    break;
	  }
	case 16 :
	  {
	    io->mod->s_opt->min_diff_lk_local = atof(optarg);
	    break;
	  }
	case 17 :
	  {
	    io->mod->s_opt->min_diff_lk_global = atof(optarg);
	    break;
	  }
	case 18 :
	  {
	    io->mod->s_opt->steph_spr = 0;
	    io->mod->s_opt->greedy    = 1;
	    break;
	  }
	case 19 :
	  {
	    io->mod->s_opt->brent_it_max = atoi(optarg);
	    break;
	  }
	case 20 :
	  {
	    io->mod->s_opt->random_input_tree = 1;
	    break;
	  }
	case 21 :
	  {
	    io->mod->s_opt->random_input_tree = 1;
	    io->mod->s_opt->n_rand_starts = atoi(optarg);
	    if(io->mod->s_opt->n_rand_starts < 1) Exit("\n. Nunmber of random starting trees must be > 0.\n\n");
	  }
	case 's':case 6:
	  {
	    if((!strcmp(optarg,"spr")) || (!strcmp(optarg,"SPR")))
	      {
		io->mod->s_opt->topo_search = SPR_MOVE;
		io->mod->s_opt->greedy      = (io->mod->s_opt->steph_spr)?(0):(1);
	      }
	    else if((!strcmp(optarg,"nni")) || (!strcmp(optarg,"NNI")))
	      {
		io->mod->s_opt->topo_search         = NNI_MOVE;
		io->mod->s_opt->random_input_tree   = 0;
	      }
	    else if((!strcmp(optarg,"best")) || (!strcmp(optarg,"BEST")))
	      {
		io->mod->s_opt->topo_search = BEST_OF_NNI_AND_SPR;
		io->mod->s_opt->greedy      = (io->mod->s_opt->steph_spr)?(0):(1);
	      }
	    break;
	  }

	case 'd':case 7:
	  {
	    if(!strcmp(optarg,"nt"))
	      {
		io->datatype      = NT;
		io->mod->ns = 4;
		io->mod->state_len     = 1;

		if(
		   (io->mod->whichmodel == LG)       ||
		   (io->mod->whichmodel == WAG)       ||
		   (io->mod->whichmodel == DAYHOFF)   ||
		   (io->mod->whichmodel == JTT)       ||
		   (io->mod->whichmodel == BLOSUM62)  ||
		   (io->mod->whichmodel == MTREV)     ||
		   (io->mod->whichmodel == RTREV)     ||
		   (io->mod->whichmodel == CPREV)     ||
		   (io->mod->whichmodel == DCMUT)     ||
		   (io->mod->whichmodel == VT)        ||
		   (io->mod->whichmodel == MTMAM)     ||
		   (io->mod->whichmodel == MTART)     ||
		   (io->mod->whichmodel == HIVW)      ||
		   (io->mod->whichmodel == HIVB)      ||
		   (io->mod->whichmodel == CUSTOMAA)
		   )
		  {
		    io->mod->whichmodel = HKY85;
		    strcpy(io->mod->modelname, "HKY85\0");
		  }
	      }
	    else if (!strcmp(optarg,"aa"))
	      {
		io->datatype              = AA;
		io->mod->state_len             = 1;
		io->mod->s_opt->opt_kappa = 0;
		io->mod->ns         = 20;
		if(
		   (io->mod->whichmodel == JC69)   ||
		   (io->mod->whichmodel == K80)    ||
		   (io->mod->whichmodel == F81)    ||
		   (io->mod->whichmodel == HKY85)  ||
		   (io->mod->whichmodel == F84)    ||
		   (io->mod->whichmodel == TN93)   ||
		   (io->mod->whichmodel == GTR)    ||
		   (io->mod->whichmodel == CUSTOM)
		   )
		  {
		    io->mod->whichmodel = LG;
		    strcpy(io->mod->modelname, "LG\0");
		  }
	      }
	    else if ((!strcmp(optarg,"generic")) || (!strcmp(optarg,"gen")))
	      {
		io->datatype = GENERIC;
	      }
	    else
	      {
		char choix;
		PhyML_Printf("\n. Unknown argument to -d option: please use `nt' for DNA or `aa' for Amino-Acids\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }

	    break;
	  }
	case 'm': case 5 :
	  {
	    int i;

	    For(i,strlen(optarg)) Uppercase(optarg+i);

	    if(!isalpha(optarg[0]))
	      {
		strcpy(io->mod->custom_mod_string,optarg);

		if(strlen(io->mod->custom_mod_string) != 6)
		  {
		    Warn_And_Exit("\n. The string should be of length 6.\n");
		  }
		else
		  {
		    Make_Custom_Model(io->mod);
		    Translate_Custom_Mod_String(io->mod);
		  }

		io->datatype              = NT;
		io->mod->whichmodel       = CUSTOM;
		strcpy(io->mod->modelname, "custom");
		io->mod->s_opt->opt_kappa = 0;
		io->mod->s_opt->opt_rr    = 1;
	      }

	    if (strcmp(optarg, "JC69") == 0)
	      {
		io->datatype              = NT;
		io->mod->whichmodel       = JC69;
	      }
	    else if(strcmp(optarg, "K80") == 0)
	      {
		io->datatype              = NT;
		io->mod->whichmodel       = K80;
	      }
	    else if(strcmp(optarg, "F81") == 0)
	      {
		io->datatype              = NT;
		io->mod->whichmodel       = F81;
	      }
	    else if (strcmp(optarg, "HKY85") == 0)
	      {
		io->datatype              = NT;
		io->mod->whichmodel       = HKY85;
	      }
	    else if(strcmp(optarg, "F84") == 0)
	      {
		io->datatype              = NT;
		io->mod->whichmodel       = F84;
	      }
	    else if (strcmp (optarg,"TN93") == 0)
	      {
		io->datatype              = NT;
		io->mod->whichmodel       = TN93;
	      }
	    else if(strcmp (optarg, "GTR") == 0)
	      {
		io->datatype              = NT;
		io->mod->whichmodel       = GTR;
	      }
	    else if(strcmp(optarg, "DAYHOFF") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = DAYHOFF;
	      }
	    else if(strcmp (optarg, "JTT") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = JTT;
	      }
	    else if(strcmp(optarg, "MTREV") == 0)
	      {
		io->datatype             = AA;
		io->mod->whichmodel      = MTREV;
	      }
	    else if(strcmp (optarg, "LG") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = LG;
	      }
	    else if(strcmp (optarg, "WAG") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = WAG;
	      }
	    else if(strcmp(optarg, "DCMUT") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = DCMUT;
	      }
	    else if(strcmp (optarg, "RTREV") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = RTREV;
	      }
	    else if(strcmp(optarg, "CPREV") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = CPREV;
	      }
	    else if(strcmp(optarg, "VT") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = VT;
	      }
	    else if(strcmp(optarg, "BLOSUM62") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = BLOSUM62;
	      }
	    else if(strcmp(optarg, "MTMAM") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = MTMAM;
	      }
	    else if (strcmp(optarg,"MTART") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = MTART;
	      }
	    else if (strcmp(optarg,"HIVW") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = HIVW;
	      }
	    else if(strcmp(optarg, "HIVB") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = HIVB;
	      }
	    else if (strcmp(optarg, "CUSTOM") == 0)
	      {
		io->datatype              = AA;
		io->mod->whichmodel       = CUSTOMAA;
	      }

	    Set_Model_Name(io->mod);

	    break;
	  }

	case 'a':case 14 :
	  {
	    use_gamma = 1;
	    if ((strcmp (optarg, "e") == 0) ||
		(strcmp (optarg, "E") == 0) ||
		(strcmp (optarg, "estimated") == 0) ||
		(strcmp (optarg, "ESTIMATED") == 0))
	      {
		io->mod->s_opt->opt_alpha     = 1;
	      }
	    else if (atof(optarg) < 1.E-10)
	      {
		char choix;
		PhyML_Printf("\n. Alpha must be > 1.E-10.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    else
	      {
		io->mod->alpha = (phydbl)atof(optarg);
		io->mod->s_opt->opt_alpha  = 0;
	      }
	    break;
	  }
	case 'b':case 10:
	  {
#ifdef ENABLE_EMPIRICAL_BAYES
	    if (atoi(optarg) < -5)
#else
		if (atoi(optarg) < -4)
#endif
	      {
			char choix;
			PhyML_Printf("\n. Branch test value must be a positive integer for bootstrap, or between -1 and -4 for aLRT branch test\n");
			PhyML_Printf("\n. Type any key to exit.\n");
			if(!scanf("%c",&choix)) Exit("\n");
			Exit("\n");
	      }
	    else
	      {
		if((int)atoi(optarg) > 0)
		  {
		    io->ratio_test       = 0;
		    io->mod->bootstrap   = (int)atoi(optarg);
		    io->print_boot_trees = 1;

		    if(io->n_data_sets > 1)
		      {
			char choix;
			PhyML_Printf("\n. Bootstrap option is not allowed with multiple data sets\n");
			PhyML_Printf("\n. Type any key to exit.\n");
			if(!scanf("%c",&choix)) Exit("\n");
			Exit("\n");
		      }
		  }
		else if (atoi(optarg)==0)
		  {
		    io->mod->bootstrap = 0;
		    io->ratio_test     = 0;
		  }
#ifdef ENABLE_EMPIRICAL_BAYES
		else if (atoi(optarg)==EMPIRICAL_BAYES)
		  {
		    io->mod->bootstrap = 0;
		    io->ratio_test     = EMPIRICAL_BAYES;
		  }
#endif
		else
		  {
		    io->mod->bootstrap = 0;
		    io->ratio_test     = -(int)atoi(optarg);
		  }
	      }
	    break;
	  }
	case 'c':case 12:
	  {
	    if ((!atoi(optarg)) || (atoi(optarg) < 0))
	      {
		char choix;
		PhyML_Printf("\n. Unknown argument to -c option: the number of rate categories must be a positive integer\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    else
	      {
		io->mod->n_catg = atoi(optarg);
		if(io->mod->n_catg < 1)
		  {
		    PhyML_Printf("\n. The number of rate categories must be a positive integer\n");
		    Exit("\n");
		  }
	      }
	    break;
	  }
#ifdef ENABLE_MULTIMAX
	case 1067:
	{
		int tmp = io->print_ml_directionality;
		int x = atoi(optarg);
		if ( (x > -1) && (x < 2) )
		{
			io->print_ml_directionality = x;
		}
		break;
	}
	case 1068:
	{
		int tmp = io->print_nni_scores;
		int x = atoi(optarg);
		if ( (x > -1) && (x < 2) )
		{
			io->print_nni_scores = x;
		}
		break;
	}
	case 1069:
	{
		int tmp = io->mod->s_opt->opt_subst_param;
		int x = atoi(optarg);
		if ( (x < 0) || (x > 1) )
		{
			io->mod->s_opt->opt_subst_param = tmp;
		}
		else
		{
			io->mod->s_opt->opt_subst_param = x;
		}
		io->mod->s_opt->opt_rr = io->mod->s_opt->opt_subst_param;
		break;
	}
	case 1071:
	{
		int tmp = io->mod->s_opt->opt_topo;
		int x = atoi(optarg);
		if ((x<0) || (x>1))
		{
			io->mod->s_opt->opt_topo = tmp;
		}
		else
		{
			io->mod->s_opt->opt_topo = x;
		}
		break;
	}
	case 1072:
	{
		int tmp = io->mod->s_opt->first_branch;
		int x = atoi(optarg);
		if (x>-1)
		{
			io->mod->s_opt->first_branch = tmp;
		}
		else
		{
			io->mod->s_opt->first_branch = x;
		}
		break;
	}
	case 1073:
	{
		int tmp = io->mod->l_max;
		double x = atof(optarg);
		if (x>0.0)
		{
			io->mod->l_max = tmp;
		}
		break;
	}
	case 1066:
	{
		io->mod->whichmodel = CUSTOM;

		Make_Custom_Model(io->mod);

		phydbl sum;

		double val1,val2,val3,val4,val5,val6;
		sscanf(optarg,"%lf,%lf,%lf,%lf,%lf,%lf",&val1,&val2,&val3,&val4,&val5,&val6);
		io->mod->rr_val[0] = (phydbl)val1;
		io->mod->rr_val[1] = (phydbl)val2;
		io->mod->rr_val[2] = (phydbl)val3;
		io->mod->rr_val[3] = (phydbl)val4;
		io->mod->rr_val[4] = (phydbl)val5;
		io->mod->rr_val[5] = (phydbl)val6;

		sum =
		  (io->mod->rr_val[0] +
		   io->mod->rr_val[1] +
		   io->mod->rr_val[2] +
		   io->mod->rr_val[3] +
		   io->mod->rr_val[4] +
		   io->mod->rr_val[5]);

		io->mod->rr_val[0] /= sum;
		io->mod->rr_val[1] /= sum;
		io->mod->rr_val[2] /= sum;
		io->mod->rr_val[3] /= sum;
		io->mod->rr_val[4] /= sum;
		io->mod->rr_val[5] /= sum;


		if(io->mod->rr_val[0] < .0 ||
		   io->mod->rr_val[1] < .0 ||
		   io->mod->rr_val[2] < .0 ||
		   io->mod->rr_val[3] < .0 ||
		   io->mod->rr_val[4] < .0 ||
		   io->mod->rr_val[5] < .0 ||
		   io->mod->rr_val[0] > 1. ||
		   io->mod->rr_val[1] > 1. ||
		   io->mod->rr_val[2] > 1. ||
		   io->mod->rr_val[3] > 1. ||
		   io->mod->rr_val[4] > 1. ||
		   io->mod->rr_val[5] > 1.)
		  {
		    Warn_And_Exit("\n. Invalid custom NT substitution rates.  All values must be greater than 0.0.\n");
		  }
		break;
	}
#endif
	case 'f':
	  {
	    if(!strcmp(optarg,"e"))
	      {
	        if (io->datatype == NT)
		  io->mod->s_opt->opt_state_freq = NO;
		else
		  io->mod->s_opt->opt_state_freq = YES;
	      }
	    else if(!strcmp(optarg,"m"))
	      {
	        if (io->datatype == NT)
		  io->mod->s_opt->opt_state_freq = YES;
		else
		  io->mod->s_opt->opt_state_freq = NO;
	      }
	    else if(!isalpha(optarg[0]))
	      {
#ifdef ENABLE_MULTIMAX
	    	if (io->datatype == NT)
	    	{
#endif
		phydbl sum;
		double val1,val2,val3,val4;

		io->mod->s_opt->opt_state_freq  = 0;
		io->mod->s_opt->user_state_freq = 1;

		/* 		sscanf(optarg,"%lf,%lf,%lf,%lf", */
		/* 		       io->mod->user_b_freq, */
		/* 		       io->mod->user_b_freq+1, */
		/* 		       io->mod->user_b_freq+2, */
		/* 		       io->mod->user_b_freq+3); */
		sscanf(optarg,"%lf,%lf,%lf,%lf",&val1,&val2,&val3,&val4);
		io->mod->user_b_freq[0] = (phydbl)val1;
		io->mod->user_b_freq[1] = (phydbl)val2;
		io->mod->user_b_freq[2] = (phydbl)val3;
		io->mod->user_b_freq[3] = (phydbl)val4;

		sum =
		  (io->mod->user_b_freq[0] +
		   io->mod->user_b_freq[1] +
		   io->mod->user_b_freq[2] +
		   io->mod->user_b_freq[3]);

		io->mod->user_b_freq[0] /= sum;
		io->mod->user_b_freq[1] /= sum;
		io->mod->user_b_freq[2] /= sum;
		io->mod->user_b_freq[3] /= sum;


		if(io->mod->user_b_freq[0] < .0 ||
		   io->mod->user_b_freq[1] < .0 ||
		   io->mod->user_b_freq[2] < .0 ||
		   io->mod->user_b_freq[3] < .0 ||
		   io->mod->user_b_freq[0] > 1. ||
		   io->mod->user_b_freq[1] > 1. ||
		   io->mod->user_b_freq[2] > 1. ||
		   io->mod->user_b_freq[3] > 1.)
		  {
		    Warn_And_Exit("\n. Invalid base frequencies.\n");
		  }
#ifdef ENABLE_MULTIMAX
	      }
	    else if (io->datatype == AA)
	    {
			phydbl sum = 0.0;
			double *vals = (double *)mCalloc(20, sizeof(double));

			io->mod->s_opt->opt_state_freq  = 0;
			io->mod->s_opt->user_state_freq = 1;

			sscanf(optarg,"%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf",&vals[0],&vals[1],&vals[2],&vals[3],&vals[4],&vals[5],&vals[6],&vals[7],&vals[8],&vals[9],&vals[10],&vals[11],&vals[12],&vals[13],&vals[14],&vals[15],&vals[16],&vals[17],&vals[18],&vals[19]);
			int qqq;
			For(qqq,20)
			{
				io->mod->user_b_freq[qqq] = (phydbl)vals[qqq];
				sum += io->mod->user_b_freq[qqq];
			}
			For(qqq,20)
			{
				io->mod->user_b_freq[qqq] /= sum;
				if (io->mod->user_b_freq[qqq] < .0 || io->mod->user_b_freq[qqq] > 1.)
				{
					Warn_And_Exit("\n. Invalid amino acid base frequencies.\n");
				}
			}
			free(vals);
	    }
#endif
	      }

	    break;
	  }

	case 'h':
	  {
	    Usage();
	    break;
	  }

	case 'i':case 9:
	  {
	    char *tmp;
	    tmp = (char *) mCalloc (T_MAX_FILE, sizeof(char));
	    if (strlen (optarg) > T_MAX_FILE -16)
	      {
		char choix;
		strcpy (tmp, "\n. The file name'");
		strcat (tmp, optarg);
		strcat (tmp, "' is too long.\n");
		PhyML_Printf("%s",tmp);
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }

	    else if (!Filexists (optarg))
	      {
		char choix;
		strcpy (tmp, "\n. The file '");
		strcat (tmp, optarg);
		strcat (tmp, "' does not exist.\n");
		PhyML_Printf("%s",tmp);
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    else
	      {
		strcpy(io->in_align_file, optarg);
		io->fp_in_align = Openfile(io->in_align_file,0);

		strcpy(io->out_tree_file,optarg);
#ifdef PHYML
		strcat(io->out_tree_file,"_phyml_tree");
#elif(MC)
		strcat(io->out_tree_file,"_mc_tree");
#elif(M4)
		strcat(io->out_tree_file,"_m4_tree");
#endif

#ifdef ENABLE_MBL
		strcpy(io->out_mixtree_file,optarg);
		strcat(io->out_mixtree_file,"_phyml_mixtrees");
#endif

		strcpy(io->out_stats_file,optarg);
#ifdef PHYML
		strcat(io->out_stats_file,"_phyml_stats");
#elif(MC)
		strcat(io->out_stats_file,"_mc_stats");
#elif(M4)
		strcat(io->out_stats_file,"_m4_stats");
#endif
	      }
	    Free (tmp);
	    break;
	  }

	case 't':case 11:
	  {
	    if ((io->mod->whichmodel != JC69) &&
		(io->mod->whichmodel != F81)  &&
		(io->mod->whichmodel != GTR))
	      {
		if ((strcmp(optarg, "e") == 0) ||
		    (strcmp(optarg, "E") == 0) ||
		    (strcmp(optarg, "estimated") == 0) ||
		    (strcmp(optarg, "ESTIMATED") == 0))
		  {
		    io->mod->kappa                 = 4.0;
		    io->mod->s_opt->opt_kappa      = 1;
		    if (io->mod->whichmodel == TN93)
		      io->mod->s_opt->opt_lambda   = 1;
		  }
		else
		  {
		    if (atof(optarg) < .0)
		      {
			char choix;
			PhyML_Printf("\n. The ts/tv ratio must be a positive number\n");
			PhyML_Printf("\n. Type any key to exit.\n");
			if(!scanf("%c",&choix)) Exit("\n");
			Exit("\n");
		      }
		    else
		      {
			io->mod->kappa = (phydbl)atof(optarg);
			io->mod->s_opt->opt_kappa  = 0;
			io->mod->s_opt->opt_lambda = 0;
		      }
		  }
	      }
	    break;
	  }
	case 'n':case 8:
	  {
	    if ((!atoi(optarg)) || (atoi(optarg) < 0))
	      {
		char choix;
		PhyML_Printf("\n. The number of alignments must be a positive integer\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    else io->n_data_sets = atoi (optarg);
	    break;
	  }
	case 'q':case 22:
	  {
	    io->interleaved = 0;
	    break;
	  }
	case 'u':case 15:
	  {
	    char *tmp;
	    tmp = (char *)mCalloc(T_MAX_FILE, sizeof(char));
	    if(strlen(optarg) > T_MAX_FILE -11)
	      {
		char choix;
		strcpy (tmp, "\n. The file name'");
		strcat (tmp, optarg);
		strcat (tmp, "' is too long.\n");
		PhyML_Printf("%s",tmp);
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    else if (! Filexists (optarg))
	      {
		char choix;
		strcpy (tmp, "\n. The file '");
		strcat (tmp, optarg);
		strcat (tmp, "' doesn't exist.\n");
		PhyML_Printf("%s",tmp);
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    else
	      {
		strcpy(io->in_tree_file, optarg);
		io->in_tree = 2;
		io->fp_in_tree = Openfile(io->in_tree_file,0);
	      }
	    Free(tmp);
	    break;
	  }

	case 'v':case 13:
	  {
	    if ((strcmp (optarg, "e") == 0) ||
		(strcmp (optarg, "E") == 0) ||
		(strcmp (optarg, "estimated") == 0) ||
		(strcmp (optarg, "ESTIMATED") == 0))
	      {
		io->mod->s_opt->opt_pinvar    = 1;
		io->mod->invar                = 1;
	      }

	    else if ((atof(optarg) < 0.0) || (atof(optarg) > 1.0))
	      {
		char choix;
		PhyML_Printf("\n. The proportion of invariable site must be a number between 0.0 and 1.0\n");
		PhyML_Printf("\n. Type any key to exit.");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    else
	      {
		io->mod->pinvar = (phydbl)atof(optarg);
		if (io->mod->pinvar > 0.0+SMALL)
		  io->mod->invar = 1;
		else
		  io->mod->invar = 0;
		io->mod->s_opt->opt_pinvar = 0;
	      }
	    break;
	  }
	case 'o':
	  {
	    if(!strcmp(optarg,"tlr"))
	      {
		io->mod->s_opt->opt_topo        = 1;
		io->mod->s_opt->opt_bl          = 1;
		io->mod->s_opt->opt_subst_param = 1;
	      }
	    else if(!strcmp(optarg,"tl"))
	      {
		io->mod->s_opt->opt_topo        = 1;
		io->mod->s_opt->opt_bl          = 1;
		io->mod->s_opt->opt_subst_param = 0;
	      }
	    else if(!strcmp(optarg,"t"))
	      {
		Warn_And_Exit("\n. You can't optimize the topology without adjusting branch lengths too...\n");
	      }
	    else if(!strcmp(optarg,"lr"))
	      {
		io->mod->s_opt->opt_topo        = 0;
		io->mod->s_opt->opt_bl          = 1;
		io->mod->s_opt->opt_subst_param = 1;
	      }
	    else if(!strcmp(optarg,"l"))
	      {
		io->mod->s_opt->opt_topo        = 0;
		io->mod->s_opt->opt_bl          = 1;
		io->mod->s_opt->opt_subst_param = 0;
	      }
	    else if(!strcmp(optarg,"r"))
	      {
		io->mod->s_opt->opt_topo        = 0;
		io->mod->s_opt->opt_bl          = 0;
		io->mod->s_opt->opt_subst_param = 1;
	      }
	    else if(!strcmp(optarg,"none") || !strcmp(optarg,"n"))
	      {
		io->mod->s_opt->opt_topo        = 0;
		io->mod->s_opt->opt_bl          = 0;
		io->mod->s_opt->opt_subst_param = 0;
	      }
	    else
	      {
		char choix;
		PhyML_Printf ("\n. The optimization parameter must be 'tlr' or 'tl' or 'lr' or 'l' or 'r' or ''.");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	      }
	    break;
	  }
	case '?':
	  {
	    char choix;
	    if (isprint (optopt))
	      PhyML_Printf ("\n. Unknown option `-%c'.\n", optopt);
	    else
	      PhyML_Printf ("\n. Unknown option character `\\x%x'.\n", optopt);
	    PhyML_Printf("\n. Type any key to exit.\n");
	    if(!scanf("%c",&choix)) Exit("\n");
	    Exit("\n");
	    break;
	  }

	default:
	  Usage();
	  break;
	}
    } // End while loop for parsing commands. . . .

  // Depricated code:
  /*   if((io->mod->whichmodel == K80) || (io->mod->whichmodel == JC69)) */
  /*     { */
  /*       if(io->mod->s_opt->opt_state_freq) */
  /* 	{ */
  /* 	  char c; */
  /* 	  PhyML_Printf("\n. WARNING: nucleotide frequencies must be set to 1/4 with this model.\n"); */
  /* 	  PhyML_Printf("\n. Type the enter key to resume the analysis.\n"); */
  /* 	  scanf("%c",&c); */
  /* 	} */
  /*       io->mod->s_opt->opt_state_freq = 0; */
  /*     } */


#ifndef PHYML
  if((open_ps_file) || (io->m4_model == YES))
    {
      strcpy(io->out_ps_file,io->in_align_file);
      strcat(io->out_ps_file, "_mc_tree.ps");
      io->fp_out_ps = Openfile(io->out_ps_file,1);
    }
#endif

#ifdef ENABLE_MBL
  if (n_l == 1) // If there's only one BL mixture, then don't bother optimizing BL mixtures.
  {	io->mod->s_opt->opt_bl_mixtures = 0;
  }

  if (n_l > 1 && io->mod->n_catg > 1)
  {
      PhyML_Printf("\n. Error: you cannot enable both the +G model and the +B model.\n");
      PhyML_Printf(". In other words, for the parameters --nclasses N --bl_mixtures B,\n");
      PhyML_Printf(". you can specify N>1 or B>1, but N and B cannot both be >1.\n\n");
      Exit("\n");
  }

#endif

  if((io->mod->s_opt->n_rand_starts)           &&
     (io->mod->s_opt->topo_search == NNI_MOVE) &&
     (io->mod->s_opt->random_input_tree))
    {
      Warn_And_Exit("\n. The random starting tree option is only compatible with SPR based search options.\n");
    }

  if ((io->datatype == NT) && (io->mod->whichmodel > 10))
    {
      char choix;
      PhyML_Printf("\n. Err: model incompatible with the data type. Please use JC69, K80, F81, HKY, F84, TN93 or GTR\n");
      PhyML_Printf("\n. Type any key to exit.\n");
      if(!scanf("%c",&choix)) Exit("\n");
      Warn_And_Exit("\n");
    }
  else if ((io->datatype == AA) && (io->mod->whichmodel < 11))
    {
      char choix;
      PhyML_Printf("\n. Err: model incompatible with the data type. Please use LG, Dayhoff, JTT, MtREV, WAG, DCMut, RtREV, CpREV, VT, Blosum62, MtMam, MtArt, HIVw or HIVb.\n");
      PhyML_Printf("\n. Type any key to exit.\n");
      if(!scanf("%c",&choix)) Exit("\n");
      Exit("\n");
    }


  if(io->mod->use_m4mod == NO)
    {
      io->mod->s_opt->opt_cov_delta      = 0;
      io->mod->s_opt->opt_cov_alpha      = 0;
      io->mod->s_opt->opt_cov_free_rates = 0;
    }

  if((io->mod->s_opt->opt_cov_free_rates) && (io->mod->s_opt->opt_cov_alpha))
    {
      io->mod->s_opt->opt_cov_free_rates = 1;
      io->mod->m4mod->use_cov_alpha      = 0;
      io->mod->m4mod->use_cov_free       = 1;
    }

  if(io->print_site_lnl)
    {
      strcpy(io->out_lk_file,io->in_align_file);
      strcat(io->out_lk_file, "_phyml_lk");
      if(io->appebr_run_ID) { strcat(io->out_lk_file,"_"); strcat(io->out_lk_file,io->run_id_string); }
      strcat(io->out_lk_file, ".txt");
      io->fp_out_lk = Openfile(io->out_lk_file,1);
    }

  if(io->print_trace)
    {
      strcpy(io->out_trace_file,io->in_align_file);
      strcat(io->out_trace_file,"_phyml_trace");
      if(io->appebr_run_ID) { strcat(io->out_trace_file,"_"); strcat(io->out_trace_file,io->run_id_string); }
      strcat(io->out_trace_file,".txt");
      io->fp_out_trace = Openfile(io->out_trace_file,1);

#ifdef ENABLE_MULTIMAX
      strcpy(io->out_annealing_file,io->in_align_file);
      strcat(io->out_annealing_file,"_phyml_sta_trace");
      if(io->appebr_run_ID) { strcat(io->out_annealing_file,"_"); strcat(io->out_annealing_file,io->run_id_string); }
      strcat(io->out_annealing_file,".txt");
      io->fp_out_annealing_file = Openfile(io->out_annealing_file,1);
#endif
    }

#ifdef ENABLE_EMPIRICAL_BAYES
  if (io->ratio_test==EMPIRICAL_BAYES && io->eb_ngens > 0)
  {
	  strcpy(io->out_eb_file,io->in_align_file);
	  strcat(io->out_eb_file,"_phyml_eb");
	  if(io->appebr_run_ID) { strcat(io->out_eb_file,"_"); strcat(io->out_eb_file,io->run_id_string); }
	  strcat(io->out_eb_file,".txt");

	  strcpy(io->out_parts_file,io->in_align_file);
	  strcat(io->out_parts_file,"_phyml_parts");
	  if(io->appebr_run_ID) { strcat(io->out_parts_file,"_"); strcat(io->out_parts_file,io->run_id_string); }
	  strcat(io->out_parts_file,".txt");
  }
#endif

  if(io->mod->s_opt->random_input_tree)
    {
      strcpy(io->out_trees_file,io->in_align_file);
      strcat(io->out_trees_file,"_phyml_rand_trees");
      if(io->appebr_run_ID) { strcat(io->out_trees_file,"_"); strcat(io->out_trees_file,io->run_id_string); }
      strcat(io->out_trees_file,".txt");
      io->fp_out_trees = Openfile(io->out_trees_file,1);
    }

  if((io->print_boot_trees) && (io->mod->bootstrap > 0))
    {
      strcpy(io->out_boot_tree_file,io->in_align_file);
      strcat(io->out_boot_tree_file,"_phyml_boot_trees");
      if(io->appebr_run_ID) { strcat(io->out_boot_tree_file,"_"); strcat(io->out_boot_tree_file,io->run_id_string); }
      strcat(io->out_boot_tree_file,".txt");
      io->fp_out_boot_tree = Openfile(io->out_boot_tree_file,1);

      strcpy(io->out_boot_stats_file,io->in_align_file);
      strcat(io->out_boot_stats_file,"_phyml_boot_stats");
      if(io->appebr_run_ID) { strcat(io->out_boot_stats_file,"_"); strcat(io->out_boot_stats_file,io->run_id_string); }
      strcat(io->out_boot_stats_file,".txt");
      io->fp_out_boot_stats = Openfile(io->out_boot_stats_file,1);
    }

  if(io->appebr_run_ID)
    {
      strcat(io->out_tree_file,"_");
      strcat(io->out_stats_file,"_");
      strcat(io->out_tree_file,io->run_id_string);
      strcat(io->out_stats_file,io->run_id_string);
#ifdef ENABLE_MBL
      strcat(io->out_mixtree_file,"_");
      strcat(io->out_mixtree_file,io->run_id_string);
#endif

    }
  strcat(io->out_tree_file,".txt");
  strcat(io->out_stats_file,".txt");
#ifdef ENABLE_MBL
  strcat(io->out_mixtree_file,".txt");
#endif


  if(io->mod->n_catg == 1) io->mod->s_opt->opt_alpha = 0;

  if(!io->mod->s_opt->opt_subst_param)
    {
      io->mod->s_opt->opt_alpha  = 0;
      io->mod->s_opt->opt_kappa  = 0;
      io->mod->s_opt->opt_lambda = 0;
      io->mod->s_opt->opt_pinvar = 0;
      io->mod->s_opt->opt_rr     = 0;
    }

  if(io->mod->whichmodel != K80 &&
     io->mod->whichmodel != HKY85 &&
     io->mod->whichmodel != F84 &&
     io->mod->whichmodel != TN93)
    {
      io->mod->s_opt->opt_kappa = 0;
    }

  io->fp_out_tree  = Openfile(io->out_tree_file,writemode);
  io->fp_out_stats = Openfile(io->out_stats_file,writemode);
#ifdef ENABLE_MBL
  io->fp_out_mixtree  = Openfile(io->out_mixtree_file,writemode);
#endif

  if(io->mod->whichmodel == GTR)
    {
      Make_Custom_Model(io->mod);
      io->mod->s_opt->opt_rr = 1;
    }

  return;
}

/*********************************************************/
